arXiv:1503.03689vl [astro-ph.GA] 12 Mar 2015 


Mon. Not. R. Astron. Soc. 000, 000-000 (0000) 


Printed 13 March 2015 (MN BTbX style file v2.2) 


The Doubloon Models 

N. W. Evans^*, J. An^, A. Bowden^ A. A. Williams^ 

^Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 OHA, UK 

^National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Road, Chaoyang District, Beijing 100012, PR China 


13 March 2015 


ABSTRACT 

A family of spherical halo models with flat circular velocity curves is presented. This includes 
models in which the rotation curve has a finite central value but declines outwards (like the 
Jaffe model). It includes models in which the rotation curve is rising in the inner parts, but 
flattens asymptotically (like the Binney model). The family encompasses models with both 
finite and singular (cuspy) density profiles. The self-consistent distribution function depending 
on binding energy E and angular momentum L is derived and the kinematical properties of 
the models discussed. These really describe the properties of the total matter (both luminous 
and dark). 

For comparison with observations, it is better to consider tracer populations of stars. 
These can be used to represent elliptical galaxies or the spheroidal components of spiral galax¬ 
ies. Accordingly, we study the properties of tracers with power-law or Einasto profiles moving 
in the doubloon potential. Under the assumption of spherical alignment, we provide a simple 
way to solve the Jeans equations for the velocity dispersions. This choice of alignment is sup¬ 
ported by observations on the stellar halo of the Milky Way. Power-law tracers have prolate 
spheroidal velocity ellipsoids everywhere. However, this is not the case for Einasto tracers, 
for which the velocity ellipsoids change from prolate to oblate spheroidal near the pole. 

Asymptotic forms of the velocity distributions close to the escape speed are also derived, 
with an eye to application to the high velocity stars in the Milky Way. Power-law tracers have 
power-law or Maxwellian velocity distributions tails, whereas Einasto tracers have super¬ 
exponential cut-offs. 

Key words: galaxies: haloes - galaxies: kinematics and dynamics - stellar dynamics - dark 
matter 


1 INTRODUCTION 

The study of spherical models is useful both for representing galax¬ 
ies and dark haloes. Even though flattening is often important, 
spherical models have provided useful insights into the behaviour 
of stellar systems. They can serve as starting points for more 
flattened models, either as initial conditions for N-Body experi- 
ments or as the lowest ord er terms in basis function expansions 
dHernouist & OstrikedI 19921) . 

An important family of spherical models discovered over the 
last twenty-five years is the d ouble-power law or y models toehnerj 
I1993I : llVemaine et^ll994h . These have a density profile that is 
cusped like at small radii, yet falls like at large radii. They 
have found ready applications in modern astronomy. They in clude 
two pa rticu larly simple and ap pealing models found earlier bv ljaffel 
( 1 19831) and lHemauis3 h99(]l) . which differ in the strength of the 
central density cusp, p ~ and respectively. All the double¬ 
power law models generate simple gravitational potentials or force- 
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laws. They also have analytically tractable distribution functions 
(DFs). This is techni cally challenging, as an integral equation ha s 
to be solved (see e.g.. lEddingtonl|l91ftlBinnev & Tremainell2008h . 
Nonetheless, DFs are useful as they encode all the properties of the 
model, enabling initial conditions to be set for N-body realisations 
of distributions of observables to be computed. 

Here, we provide another very simple family of spherical 
double-power law halo models - the doubloon models. They are 
motivated by the flatness of galaxy rotation curves, and so the mod¬ 
els have a density that falls like either in the inner parts or the 
outer parts. This generates a regime in which the rotation curve is 
roughly flattish. Explicitly, we consider the potential-density pair 



(l + p)aP + U’ 
^ AnGr^^e (aP + rP)^ 


( 1 ) 


Then, the rotation curve of the model is 


Vcir) = V| 


('Ar'" 


( 2 ) 


Here, V is the amplitude of the flat rotation curve, whilst a is a 
scalelength. The density is positive everywhere provided p > -1, 
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whereas it is monotonically decreasing (and hence astrophysically 
realistic) provided p < 2. When p = 0, the model degenerates 
into the singular isothermal sphere (i.e. p oc but the strict limit 
is scaled to Vc = V/ V2), known for over a century thanks to th e 
labours of J. H. Lane and R. Emden (see also lChandrasekhaJ 1 9391) . 

The properties of the family divide neatly into two, as shown 
by their central and asymptotic behaviour. If p > 0, the density falls 
off like as r —> oo and the rotation curve is flat asymptotically 


Vc(r) =^Vx 


{r/ayl^ 


(r <K a) 
(r » a)’ 


(3) 


which we shall subsequently refer to the outer branch. The poten¬ 
tial of the outer branch decreases from = 0 to if/(oo) = -co as 
r increases, and the models possess the central density cusp like 
(except for p = 2). Some of the usual suspects are rep¬ 
resented in the family. When p = 2, the model is the spherical 
limit of Binney’s logarithm ic potential llBinne van 1 19931: 

iSinnev & Tremainal2008h . When p = 1, the model has a 1/ r cusp 
and has recently been discussed bv lEvans & Williams! (|2014|) . 

The rotation curve for models with p < 0 tends to a finite 
value at the centre 


Vc(r)=^Vx 


(fl/r)l'’l/2 


(r <K a) 
(r » a)’ 


(4) 


Vc/V 




but falls off like as r —> oo. Henceforth these models will be 
referred to as the inner branch. The density always has an isother¬ 
mal cusp (~ r^^) at the centre and decays like as r —> oo 

(with the exception of the case p = -1), whereas the potential runs 
from i/^(0) = oo to i/f(oo) = 0. The model with p = -1 and was 
introduced bv ijaffel (Il983h as a representation of elliptical galaxies 
and has finite mass (M = aV^IG). The remaining members of the 
family have rotation curves which fall off less steeply than the Jaffe 
model. 

Both branches of the doubloon family are useful. Eor exam¬ 
ple, models on the inner branch are helpful in studies of the outer 
parts of galaxies, when the flat rotation curve gives out and the den¬ 
sity of the dark matter begins to fade. Evidence, for example, from 
the Sagittarius stream suggests that the rotation curv e of the Milky 
Way is flat out to 50 kpc, and then begins to fall dOibbons et all 
l2014h . Models on the outer branch are useful for studying the in¬ 
ner parts of galaxies, in the regime of the flat rotation curve. Their 
central cusps make them desirable models of dark matter haloes, in 
acco rd with predictions of dissipationles s theories of galaxy forma¬ 
tion dMo. van den Bosch & Whitj|201(]|) . 

The paper is arranged as follows. Section 2 studies the self- 
consistent model, and gives families for the distribution functions 
for the total (dark and luminous) matter. These are self-consistent 
models and so are useful both for setting up the initial conditions 
for N-body experiments and for studying the kinematic properties 
of the dark and luminous matter. The remainder of the paper stud¬ 
ies tracer populations moving in the doubloon potential. The trac¬ 
ers might represent stellar populations in elliptical galaxies or the 
spheroidal components of spiral galaxies. In particular. Section 3 
looks at the Jeans solutions of tracer populations (with power-law 
or Einasto profiles), whilst Section 4 studies the distributions of 
high velocity stars. Both applications are motivated by the data on 
halo stars in the Milky Way, which has increased in quality and 
qua ntity in recent years thanks t o surveys like SPSS and RAVE (se e 

e.g.. lSmith et al.l200ll200^fl!Bond et al.l2010l:|Piffl et al.l20l4 . 


Figure 1. Rotation curves Vc (upper) and isotropic velocity dispersion CTp 
profiles (lower) for the models p = 2 (black), p = 1 (red), p = 1/2 (yellow), 
p = -1/2 (green) and p = — 1 (blue). 


2 THE SELF CONSISTENT MODEL 

The phase space distribution function ( DF) may depe nd on the in¬ 
tegrals of motion only, as first realised bv ijean 3 J1919I) . Eor a spher¬ 
ical potential, the integrals may be taken as the binding energy per 
unit mass E and the modulus of the angular momentum L, namely 

E = -{{vl + vl+vl) + ilr, L^ = r\vl+vl). (5) 


Isotropic DEs depend only on E, whereas anisotropic DEs depend 
on L as well. Here we look for constant anisotropy DEs; that is, 
the anisotropy parameter /J = 1 - {v\)l{vl) takes a constant value. 
The parameter p may take values in the range -oo < /? < 1. 
When yS —> 1, the model is built from radial orbits, whilst when 
P -oo, it is made from circular orbits. Constant anisotropy DEs 
are widely-used because of their simplicity. However, it is impor¬ 
tant to acknowledge that there is no underlying physical justifica¬ 
tion. Simulations of the growth of galaxies suggest that DEs are 
typically isotropic in the center (P x 0), but become more radiall y 
anisotropic (p x 1/2) in the outer parts toansen & Moor j|200d) . 
It is nonetheless reasonable to expect that constant anisotropy DEs 
provide good approximations in certain regimes, such as the central 
parts or the outer periphery. 

The velocity dispersions of the constant anisotropy models are 
found by integrating the Jeans equation, namely 




( 6 ) 


and (Ug) = (1 -j8)(uJ). For the self-consistent doubloon models in 


© 0000 RAS, MNRAS 000, 000-000 
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equation Q on the outer branch, assuming 0 < p < 2(1 -/3), 


2.1 Outer Branch 






(\ + p)aP + rP 


d^'->^-P)(aP + rPf 






(7) 


where A = 2p“'(l -j3)- I >0 and ij = gPHgP + rP), whilst B ,(a, b) 
is the incomplete beta function jAbramowitz & Stegunll96^ § 6.6; 
lOlver et ^1201(1 § 8.17). For all members, the velocity dispersion 
tends asymptotically to ^ Vy-^2(1 -P) as r ^ co. As r 

0 , it behaves like (vf.) ~ piim(p, 2 - 2 /i-p)^ exception of the 

case p = I - p when (v^) ~ rP log r. The radial velocity dispersion 
velocity tends to zero at the center, except for when p = \ — pjl. 

For the inner branch models (-1 < p < 0), we find 






\+{l-\p\)x 


0^ 


2{1-I3) 


(1 + x)2^B 1 (2 + ^,-^)-- 


( 8 ) 


where ^ = 2|pr‘(l -^6) = -1 - d > 0 and x = (a/r)" = (r/^. 
This is finite at the centre forjS < 1, as limr^o(i'r)'^^ - F/-^2(l - /?), 
whilst it decays like {v^) (V^/2)(l - P + |p|)“*(a/r)''’' as r —» oo. 

Plots of the isotropic (p = 0) velocity dispersion are shown in Fig[T] 
as an illustration of some of these properties. 

The key to the simplicity of the doubloon models is that if)(r) 
can be easily inverted, namely 


A physical model must have an everywhere positive DF. Since p ~ 
d^P as r ^ 0 for the self-consiste nt model with p > 0 , the cusp 
slope-central anisotropy theorem of lAn & Evans! ( l2006all indicates 
that the constant anisotropy DF is physical only if 2/1 < 2 - p. The 
constant anisotropy model for p > 0 has the DF expressible as the 
sum over the exponentials of p£ where £ = E/V^, namely 


_ (l-/l)3/2-d a2,-2y2,-, 

47t5/2F(1 -P) GL^ 




(d-l); (1-Hp)(d) 


0 - 1 ) 




pj>6 


(13) 


where (A)j = O/po ('1 + 01® the Pochhammer symbol and F(x) is the 
gamma function. Also note d = 2(1 - P)Ip - 1 > 0. For p = 1 (the 
halo model wi th the 1/r density cusp), t his reduces to the expres¬ 
sion given by lEvans & William^ j2014h . Experimentation shows 
that the sum over the exponentials converges rapidly, so this is a 
practical way to compute the DF 

There is a particularly simple case that is worthy of note. If 
d = 0 (i.e. p = I- p/2), then the DF reduces to just a sum over two 
isothermals, multiplied by a power of the angular momentum 


F = 


V'-PLP-^ 
47t5/2F(p/2) GaP 



ePS 



(14) 


r(i/r) = fl|exp(-||^j-l| . 

This means that plip) can also be easily constructed 


p(d) = 


q2PIv2 pg(p+2)^/v2 
47tGa2 _ gpiA/v2^2/p-i 


(9) 


( 10 ) 


From this, we can use Eddington’s (1 916) formula to derive th e 
isotropic DF via Abel transforms (see lBinnev & Tremainell2008h . 
For anisotropic DFs, we similarly construct the augmented density 


p(<A) = (2rYp(iA) = 


24-2 y2 


gPiA/v2 p^2pil//V^ 


nGa^^^ 4) ^g-p((,/v2 _ j^2(i p)ip 1 
Then the constant anisotropy DF has the form of: 


F(E,L) = 


fE{E)Hp{L^) 
(27T)3/2 ’ 


9 I- (j6<l) 

Hp{L^) = F(1 -P) ^ \ 

W) ip=i) 


( 11 ) 


( 12 ) 


where 6{x) is the Dirac delta function, whilst p(i/r) at a fixed p 
is given as an integral transformation of fsiE). The explicit rela- 
tionship between pf i p) and fsiE) is f ound in the l i teratu re; e.g., 
IWilkinson & Evansl j 19991 eq. 2) or lEvans & At] j2006L eqs. 2 
& 3). We outline a general algorithm to find /gfE) given g(tf/) in 
Appendix lAl which is an elaboration of Eddington’s inversion for 
isotropic DFs in terms of Abel transforms. 

We note that the outer branch gives simpler DFs than the inner 
branch. For the outer branch, the DF is always a sum over isother¬ 
mal or Maxwell-Boltzmann distributions. Usually, the sum is in¬ 
finite, but, for some special cases, the sum is finite. For the inner 
branch, the DF is a sum over inco mplete gamma functions (cf . 
lErdelvi et alJll95^ vol. 2, chap. IX; lAbramowitz & Stegun|[T964 


§ 6.5; Olver et aDbOIoL chap. 8). Under some circumstances, the 
sum is finite and over Dawson’s integrals (which are equivalent to 
error functions). Although we give the general solutions, we point 
out some of the simple cases along the way. 


These are amongst the simplest DFs known, and some particular 
cases of this family have already appeared in the literature. So, for 
Binney’s logarithmic model (p = 2), 


F = 


1 




iexp(^) + V2exp(^) 


(15) 


which re duces t o the isotropic DF comprised of two isothermals 
given in lEvansI ( Il993h . For the halo model with the 1 jr density 
cusp (p = 1), it reduce s to the simple anisotropic DF found in 
lEvans & William j(l2014h . 


F = 


1 

n It 

1 /2£\1 

TxAGaF 


2)^4“Kv?)| 


(16) 


For all p, the velocity dispersions corresponding to the DF of equa¬ 
tion d are analytic and given by: 


{v]) = 


(2 + p)aP + 2rP 
2 p {I + p)aP + rP 


, {2 +p)aP+2rP 

(t'fl) =-. 

4 (1-1- p) aP + rP 


(17) 


This provides a simple solution to the Jeans equations for the ve¬ 
locity dispersions. 

We remark that (i) the same simplicity occurs for the hyper- 
virial model^ when the cusp slo pe p and anisotropy p are related 
P = 1 - p/2 I Evans & Anll200^ and that (ii) this combination of 
cusp slope and anisotropy is at the extreme permitted by the cusp 
slope-anisotropy theorem. 

If /I = 1 (i.e. p = I - p) on the other hand. 


F = 


pl/2+pyl-2p^2(p-l) 

47t5/2F(p) Gd^P 


(e2'’fi + (I-(-p)£(|)'"'’e2>s). (18) 

3=3 


^ Hypervirial models satisfy the virial theorem locally. The most celebrated 
example is the IPlumme model , for which the property of hyper- 

viriality was established by lEddingtorJ jl916h . lEvans & ArJ j2QQ ^ foun d 
a family of models with this property - See also lAn & Evansl f2Q06bl). 
The pr oper ty of hypervirialit y has been studied theoretically bv llguchi et all 
ilOOd) and lSota et al.l i2008h 
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with the velocity dispersions expressible analytically to be 
V^rP 




(I + p)aP + rP 


\+p 


(l + —)^og(l + —— 
\ aP / \ rP / aP 


2 +'ip 

2p 


(19) 


Finally, we give the isotropic (fi = 0 and so d = 2/p - 1) or 
ergodic DF that depends only on energy, which is, for 0 < p <2 


F = 


1 


1 + 


expl 


(w) 




, (l+p)(2-l),._. 

- + 


1=1 








JpS 


( 20 ) 


For p = 2, this simplifies to equation 03. whilst if p = 1, this 
reduces to 


F = 


1 


An^FGa^V 


exp(^) + (l+p)|;(^)'exp(^) 

1=3 


( 21 ) 


The isotropic velocity dispersion resulting from the ergodic DF is 




V^r" 


(\ + p)aP + rP 


r^V-p)^aP + rPf 2 


2/2 2 \ 1 


where y = a'V(o^ + r’’), with the particular cases of 


(ii?) 




2a^ + E 
2 ia^ + F- 
W \2{a + rf 


(P = 2), 


f2(a + r)^ / 5a + 4r^ 


2 a + rl 


( 22 ) 


(23) 


As r —> 0, the velocity dispersion tends to zero (unless p = 2, for 
which (vj) V^li), whilst it tends to (uj) 1/^/2 as r ^ oo. 


which therefore provides a simple radially anisotropic DF for all 
the doubloon models on the inner branch. For p = -1, 

exp(£) - 3 exp(-£) + 2 exp(-2£) 

^ -’ ^25) 

which is the DF for the Jaffe sphere. This should be particularly 
useful in setting up initial conditions for N-body experiments. 

All the models on the inner branch exhibit an isothermal cusp 
and so it is technically possible to set up the model composed en¬ 
tirely of radial orbits, although the resul ting models will in gen¬ 
eral be prey to the radial orbit ins tability l lFridman & Polvachenkd 
ll984[Mrner & Papaloizoull987h . The DFs of the radial orbit mod¬ 
els are given in the closed form; 

F=^^y\f\(l + \p\)Dp{y[\i^)-y/2\p\iDp[yl^)\, (26) 

utilizing Dawson’s integrals Isee lBinnev & Tremainell200^ . though 
here we use D± instead of F± to avoid confusion with the DF): 


DJ^x) = exp(-HX^) 



exp(±t^) df. 


(27) 


This is actually closely related to the error function; that is, D+(x) = 
(V^/2)e“'’^erfi(x) and D_(x) = (V^/2)e'’^erf(x), where erfi(x) = 
i“*erf(ix) is the imaginary error function. The radial velocity dis¬ 
persions corresponding to the pure radial orbit models are express¬ 
ible using only elementary functions 

. V^ (l-^x)^log(l+x-‘)-x-l-|p|/2 

Ipl l + (l-|p|)x ’ ^ ’ 


where x = (alr)P = {rlajp\ which behaves like (vj) log(a/r) 
as r —> 0 and (vf) V^(2|p|)“'(a/r)''’' as r —> oo, respectively. 


2.2 Inner Branch 

The inner branches possesses more complicated DFs t han the outer 
branch. This may be guessed from the properties of the |jafrel ( ll983l) 
model, wh ose isotropic DF, first foun d by Jaffe and subsequently 
reported in iBinnev & Tremainel (l2008h . is already a sum of special 
functions (viz. Dawson’s integral). 

Once we expand g(^) in a power-series in e,Fp'*l'P~ < 1, the 
constant anisotropy DF in general is expressible as a sum over 
the incomplete gamma functions. The sum terminates after a fi¬ 
nite number of terms if 2|p|“‘(l — fi) is a non-negative integer. If 
j8 is a half-integer, the actual operation reduces to ordinary deriva¬ 
tives (rather than Abel transforms or fractional derivatives) and so 
the result is eventually expressible by means of a rational function 
of exponentials. If /? is an integer, the resulting incomplete gamma 
functions are reducible to the error function (or equivalently Daw¬ 
son’s integral). If \p\ = 2(1 -j8)/n < 1 where « is a positive integer, 
the final expression is resolved into a finite sum over such func¬ 
tions, as is the case for the isotropic Jaffe model. 


2.2.1 Models with radially biased orbit distributions 

The simplest constant anisotropy DFs for the inner branch models 
are obtained for 13= 1/2. These models are of widespread phys¬ 
ical applicability, as numerical simulations suggest that these are 
characteristic of dark m atter haloes, at least in the outer parts (e.g., 
iHansen & Moorell2006h . Using the notation £ = E/V^, we find: 

(glpie - , 

F = ^ -(1 -H Ippe^'"'® - 2|ppe-2''’'®), (24) 

(27t)3GaL V r 


2.2.2 The Jaffe model 


In fact, the Jaffe model deserves particular attention, as it is of- 
ten used a model of a dark halo or a n elliptical galaxy (see e.g., 
lKochanel3ll996l : lOerhard et al.lll998h . The radial velocity disper¬ 
sion of the constant anisotropy Jaffe model is given by 


ivl) = V- 


2 r^F + ,-)2 


a^F-H) 


B^ (5-2/1,2/3-2), 


(29a) 


which is reducible to an elementary function if 2/3 is an integer - in 
particular, if n = 2 - 2/3 is a non-negative integer. 


(V^r) 

V2 


®„(x) = £ 


(l+n)(2 + n) 
2 

" (- 1 )*-' 


122 ^ (n + fe)x* 


l-v[log(l + i) + 2 


3 + « 2 + n r 


2 a 

(-If 

krJ 


(29b) 


An analytic DF of the Jaffe model with /I = 1/2 is already 
provided in equation ( 125b . Similar analytic DFs of the Jaffe model 
may also be obtained for all half-integer values of the anisotropy 
parameter. For instance, if/3 = -1/2, 


(e®- 1)^(9+ 7e-®+4e-2e)L 
(27t) 3 Ga31/2 


(30) 


For an integer value of the anisotropy parameter, the constant 
anisotropy DF of the Jaffe model is expressible as a finite sum over 
Dawson’s integrals. In particular, the DF with only radial orbits is 
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5 


whilst the Jalfe models with -/i = n 6 jO, 1,...): 


j=-2 


2«+5/27t3„!G 

Q-2(l+n) y-(l+2n)jy.n r^i^+l) 


,11/2 


7 + 2 


D 




2n+5/27^5/2„|G 


z 

J=i 


{-iyj"+2 (211 + 4 

7 + 2 


- (-!)"(« + 2)e-^erfi(V5) + {-2fyl2&-'^e.Yfi{\l2£) 


(32) 


where D(, =Q and is the binomial coefficient. The isotropic JafFe 
model is included as the special case {n = 0) 

^ Z7+(V2£)-V2Z)+(V£) + D_(V2£)-V2Z)_(V£) 

which was first giv en EiS l ll983h and is also repeated in 
[Binnev & Tremaind f20()it . 

The JafFe models given here all have co nstant anisotrop y. They 
may be contrasted with the models found bv lMerritl ( Il985ah . which 
have isotropic centers and strongly radially anisotropic (,0 —> 1) 
outer par t s. The se are derived using the inve r sion int roduced by 
lOsipkovI h979l) and popularised by iMerrid ( Il985bh . Although 
the transition from isotropy to radial anisotropy is desirable, the 
Osipkov-Merritt models unfortunately provide rather too extreme 
radial anisotropy in the outer parts. 


3 FLATTENED TRACER POPULATIONS: JEANS 
EQUATIONS 

For applications in which dark haloes are represented as doubloon 
models, we are primarily interested in the properties of flattened 
tracer populations of stars, whose kinematics are accessible to ob¬ 
servation. The flattening is usually described by assuming a density 
law stratified on similar concentric spheroids with an axis ratio q. 
Tracers in haloes are often modelled with power-laws or broken 
power-laws (e.g.. IWatkins et aill2009l : [Peason et al.ll201 ih . which 
have asymptotic behaviour 

Ppi(T 6) — AmT^, rri^ - r“(sin^ 6 + cos^ 6). (34) 

Another commonly used profile for luminous tracers, whether in 
stell ar haloes or elliptical galaxies, is that due to J. Einasto (see 
e.g.. lEinasto & Haudlll989riDeason et alj201lh. though i t has also 
been used for the dark matter distribution Joraham et al.ll2006h . It 
has asymptotic behaviour 


PEin(T 0) - A exp(-im‘/"), (35) 

where i and n are constants. This is the exponential profile if « = 1 
and the de Vaucouleurs-like profile if n = 4. A convenient general- 
izat ion that incorpora tes both families is the hyper-Einasto profile 
(cf. lAn & Zh^l2013h 

PHyEm(T &) - AivT^ exp(-im‘/"). (36) 

We shall solve the Jeans equations for both flattened power-law and 
hyper-Einasto tracers in dark matter haloes described by doubloon 
potentials. In an axisymmetric model, {v^v^) = {vgVip) = 0 by the 
symmetries of the individual orbits, which then leaves two non¬ 
trivial Jeans equations for the tracer density p,(r, 6) in terms of the 
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spherical potential i/'(r): 

9pt<i^?) 


9r 


1 ^pAvrVe) A 
r do r 


+ l^[2{vl)-{vl)-{vl) + {v,ve)co\B) 


dif/ 




— Pi “ “Pt 7 ’ 
or a’’ -I- rP 


dpt(UrV0} 1 ^PtiVg) fh 

dr r do r 


+ — [3<uas> + (.dl) - <!/^»COt0] = 0. 


(37) 

(38) 


These two relations between the four stresses p{v(.), p{vrVg), piv^), 
p{vp must be satisfied at any point in the model. 

There is now increasing evidence from the stellar halo of our 
Galaxy that t he velocity e l lipsoid is spherically aligned. This was 
first noted by ISmith et^ ( l2009al) . who studied the kinematics of 
halo subdwarfs in the Sloan Digital Sky Survey Stripe 82, for which 
there is multi-epoch and multi-band photometry per mitting the 
measu rement of accurate proper motions. Subseauentlv. lBond et al.l 
l l2010h used ~ 53000 halo stars with r band magnitude brighter than 
20 and proper motion measurements derived from Sloan Digital 
Sky Survey and Palomar Observatory Sky Survey astrometry to ex¬ 
tend this result over a quarter of the sky at high latitudes. Although 
the tilt of the velocity ellipsoid in elliptical galaxies is not known, 
galaxy mo delling sugges ts that it is aligned on spheroidal coordi¬ 
nates (e.g.. lBinnevll2014l) . which become spherical at large radii. 
In other words, there is much motivation for investigating spherical 
alignment, (v^vg) = 0, as it holds good for the Milky Way’s stellar 
halo and for the outer parts of elliptical galaxies. 

Additionally, there is good evidence from the ki nematics of 
stars in the Milky Way’s stellar halo (Vg) x (vp. ISmith et'H] 
(l2009bl) fou nd = 82 ± 2 km s“‘ and = 11 ±2 km s“*. 
iBond et ^ ( l201oli claimed that the semi-axes are invariant over the 
volume probed by their much larger sample and found = 

85 ± 5 km s“' and = 75 ± 5 km s“'. It is interesting that 

the two angular semi-axes are almost the same. If the tracer popu¬ 
lation has a DE depending on E and L only, then (ug) = (u^). We 
accordingly make this assumption, which has good theoretical and 
observational motivation, so that the Jeans equations become 


dr 


= rp^{2{vl) - vl). 


9(Pt<fe» 

do 


(39) 


We see that this set of assumptions has closed the Jeans equations, 
which may now be integrated with suitable boundary conditions. 
Integrating the angular equation, we find pt(i'e) being an arbitrary 
function of r. As the boundary condition, we next assume that the 
velocity dispersion on the equatorial plane (0 = n/2) is a constant 
fraction K of the rotation curve, that is. 


(dg) = {vp = Kv((r) 


Pijr, 71/2) 

Pt(r, 0) 


(40) 


Assuming that p{v(:) —» 0 as r —» oo, we then have the full axisym¬ 
metric solution of the Jeans equations as 


(d?) = G 0) X”7t/2)] (41) 

This gives a one-parameter family of solutions of the Jeans equa¬ 
tions for axisymmetric densities with spherically aligned velocity 
dispersion tensors. Algorithms for cy lindrically al igned Jeans so- 
lutions are known dCappellaril l2008l) . although as iBinnevI ( l2014h 
points out they are somewhat contrived. In cosmogonies in which 
galaxies are built from hierarchical merging, stellar material on 
nearly radial orbits fell in to deepening dark matter potential wells, 
and so spherically aligned Jeans solutions are much more natu¬ 
ral. A general, though rather complicated, algorithm for spheri¬ 
cally aligned Jeans solutions for flattened potential-density pairs 
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Figure 2. Velocity ellipsoids for flattened (q = 0.6) Einasto profiles (upper panels) and power-law profiles (lower panels). The left column is for doubloon 
halos with p = I (outer branch), the right column for p = —I (inner branch). All the doubloon models have V = 200 km s“* and a = 50 kpc. The Jeans 
solutions correspond to the case K = 1/40. Notice that the flattened Einasto profiles do not have approximately invariant shapes, as the long axis of the velocity 
ellipsoid changes from radial to azimuthal at high latitudes. 


is known llBacon et alJll983l : lBaconlll985l) . The assumption of a 
spherical potential, though, makes our algorithm much simpler for 
flattened tracer densities. 

A necessary condition for everywhere positive stresses is that 
0 < K < 1/2. In practice, the range of physical K values is 
much more constrained, though established easily enough by nu¬ 
merical integration. Fig.|2]shows velocity ellipsoids for power-law 
and Einasto prohles representing the stellar halo of the Milky Way. 
The Einasto profile has « = 1.7 and an effective radius of 20 kpc, 
the power law profile falls with d = -4 beyond a core radius of 
0.6 kpc. Both ar e inspired by fits to the st ellar halo of the Milky 
Way d iscussed in lPeason et alj ( l201 lLl2014h and lEvans & WilliaiM 
( 120141) . The l eft column shows each tracer in a doubloon model 
with p = -1 ) |jaffelll983h . the right with p = 1 dEvans & Williams! 
l2014h . It is interesting to observe that the shape of the velocity ellip¬ 
soids is primarily controlled by the tracer density, with the power- 
law or Einasto profiles each generating similar Jeans solutions in 
different doubloon potentials. However, the shape of the velocity 
ellipsoids for power-laws tracers always has the radial velocity dis¬ 
persion exceeding the angular velocity dispersions, so that the ve¬ 
locity ellipsoids are always prolate spheroids. This is not the case 
for the Einasto profiles, in which the azimuthal velocity dispersions 
exceeds the radial on approach to the poles (6 = 0), and so changes 
from prolate to oblate spheroidal in shape. 

Although IPeason et alJ (l20lT|, l2014h found either power-law 


or Einasto profiles equally good fits to the starcount data, it is obvi¬ 
ous that the kinematics provides a powerful discrimin ant. The fact 
that th e velocity ellipsoid shape is spherically aligned dSmith et'H] 
l2009ah and (to first o rder) shape invari ant over the Sloan Pigital 
Sky Survey footprint dBond et al.ll2010l) seems to rule out Einasto 
profiles. We plan to return to detailed Jeans solution fits to the kine¬ 
matics of the stellar halo in a later publication. 


4 TRACER POPULATIONS: DISTRIBUTIONS OF HIGH 
VELOCITY STARS 

The high velocity stars of the Milky Way are distinct fro m the 
hyper-velocity stars. The central black h ole may eject stars dHillsI 
Il988l : l^ & Tremainel2003l : lLevirdl20()^. which M e often unbound 
and moving on highly radial orbits i Brown et alj|2014l) . These are 
the hyper-velocity stars, which are a separate population and do not 
form part of the steady-state stellar halo. 

By contrast, the high velocity stars are the highest energy, but 
bound, members of the halo. The form of the distribution func¬ 
tion at the highest energies is accessible to observational scrutiny 
and can in principle provide information on the behaviour of the 
potential at the edge. For example, in the Milky Way, the distribu¬ 
tion of high-velocity stars fro m the halo is alread y available locally 
thanks to the RAVE survey dSmith et ^ 1 20071 : IPiffl et al] l2014l : 
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Figure 3. The probability density function of radial velocities near the escape speed for power-law (left) and Einasto (right) tracer densities. The logarithm 
of the probability density is plotted against 1 - p/yesc- The red line is for an isotropic power-law tracer with (5 = 5 in a doubloon potential with p = 1 (outer 
branch, the Evans & Williams (2014) model). The black and blue lines show tracers in a doubloon potential with p = -\ (inner branch, Jalfe model) with 
^ = 5 and 4. Full and dashed lines show isotropic (J3 = 0) and radially anisotropic (J3 = 0.5) models. Although all four models are shown in the right panel, the 
super-exponential cut-olf makes them virtually indistinguishable. 


iHawkins et al.|[2QT^. One of the earlies t investigations into the es¬ 
cape speed, Leonard & Tremain a l ll990h introduced the simple and 
attractive ansatz that 


\f(v) oc - I')'^ (v < t)esc) , .-x 

\/(i;) = 0 (V> t)e.c), 

where f(v) is the distribution of space velocities near the escape 
speed fesc and C is a constant. This ansatz, which is exact for 
isotropic power-law DFs , has held up su rprisingly well over the 
last quarter of a century dPiffl et al.ll2014l) . However, the next few 
years will see the Gaia-ESO and LAMOST surveys, as well as the 
Gaia satellite, substantially improve our knowledge of the distribu¬ 
tion of high velocity stars as a function of distance within 20 kpc 
of the Sun. Sample sizes of hundreds or even thousands of high 
velocity stars will become available, and deviations from equation 
(ES can be probed. Accordingly, we proceed to derive the form of 
the velocity distribution at the highest energies for power-law and 
hyper-Einasto tracers (defined in ea.l36t. 


4.1 Outer Branch 


The highest velocity stars in the outer branch models correspond to 
the limit E —> -oo, and we find r(E) —» oo in the same 

limit. If the tracer density asymptotically becomes a power-law like 
p as r —» oo, then for E —» -oo. 


g - exp(^L|^). 


(43) 


We find that the asymptotic form of the constant anisotropy DE is 


2 ^A . 

(S-2p\ 


f (S-2P)\E\ 

(27t)3/2a«-2/jl 

[ V2 J 

' F(1 -/3) 

[ V2 


(44) 


for E —» -oo. So, an isotropic power-law tracer (Ji = 0) has an 
isothermal or a Maxwellian distribution of velocities. Even in the 
presence of anisotropy, the distribution of radial velocities remains 
a Maxwellian. The red line in the left panel of Fig. [3] shows the 
probability density function (PDF) of radial velocities near the es¬ 
cape speed derived from eq. El for the model with d = 5 and 

j8 = 0. 

If the tracer population has an hyper-Einasto profile, then g == 


2^Ar^ 3 exp(-.?r'^") as r ^ oo and the DF asymptotically becomes 


2f‘A 

(Ja'7' 

L-^ 

(27t)3/2a«-2/J 

1 nV^ 

X exp 

r(i -p) 

1 _EL / 

-jfl" e»r2 —\s — 2p 


3-2^ \ |£| 
2n I 


(45) 


as £ —> -oo. The distribution of space velocities is no longer 
Maxwellian, but rather is a Maxwellian modulated by a super¬ 
exponential cut-off. 


4.2 Inner Branch 


The forms of the high energy tail of the velocity distribution 
change if the dark matter density falls off like or faster, as 
we now show. For p < 0, Taylor expansion now shows that 
r{E) ^ —» oo in the limit £ —> 0. For power- 

law tracers, then as r —> oo and £ —> 0, 


g - 2l^Ar^-^ - 2l^Aa^H ^ f (where ^ ) (46) 

\ W / \p\ 

which is also a power-law and so the asymptotic form of the DF is 


2^A FQr+l) £-^^ 

(27t)3/2a'5-2/s T{x-il2+P)T{\-li) \ I 


(47) 


for £ —> 0. In other words, the space velocity distribution of a 
power-law tracer falls asymptotically like a power-law at the high¬ 
est energies. The black and blue lines in the left panel of Fig. |3] 
show the PDF of radial velocities derived from equation I47t for 
the model with 6 = 5 and 6 = 4. Full lines are isotropic = 0), 
dotted lines radially anisotropic = 1/2). 

For an Einasto tracer, a similar calculation yields 


£ ^ 


2l^A L-^ 
(27t)3/2a«-2/i T{\-P) 


^ , 3 o 5-2B 3,-28 

n|p|£/ / 



(48) 


as £ —» 0. So, the distribution of space velocities falls like a power- 
law with a super-exponential cut-off. The same four models (6 = 5 
and 4, and p = 0 and 1/2) are shown in the right panel of Eig. [3] 
The super-exponential cut-off causes all four model to lie on top of 
each other. 
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5 CONCLUSIONS 


We have presented details of a new family of spherical models, 
which have properties suitable for mimicking galaxies with flat- 
tish rotation curves. The models may have a rotation curve which 
attains a finite value at the centre and falls on moving outwards. 
The archetype is the Ijaffd ( Il983h model. Alternatively, the mod¬ 
els can have a rising rotation curve in the inner parts, which flat¬ 
tens asymptotically. Here, the archetype is the spher ical logarith¬ 
mic model popularised bv [Binnev & Tremainel ( l2008h . The family 
also includes the singular isothermal sphere with an p ~ den- 
sit y cusp at the center, as wel l as the halo model recently discovered 
bv lEvans & Williams! ( 12014l) which has an cusp. In general, the 
family includes both cored and cusped members, and so can rep- 
resent the range o f cusps found in numerical simulations (see e.g., 
iMoore et Ii]| 19991) . 


The halo models presented here are spherical. Cosmologi¬ 
cal simulations suggest that dark halos are typically flattened and 
mildly oblate, with a ratio of long axis to short axis of p = 0.8- 
0.9 (see e.e.. lAbadi et al. IboiOl : Ipeason et ^1201 ll : IZemp et^ 
l2012h . For the Milky Way Galaxy, there are several lines of ev¬ 
idence suggesting that the dark halo may be nearly spherical. 
First , the fits to the tidal stream GD-1 usin g different meth¬ 
ods dKoposov et alj|20l3 : ISowden et al.ll2014h show that the to¬ 
tal Galactic potential (disk and halo) at radii 15 kpc is consis¬ 
tent with modest flattening. Secondly, the kinematics of halo sub¬ 
dwarfs in the Sloan Digital Sky Survey Stripe 82, for which there 
is multi-epoch and multi-band photometry permitting the measure¬ 
ment of accur ate proper motions, is consistent with a nearly spher¬ 
ical potential dSmith et al.ll2009i) . Thirdly, the kinematics of the 
Sagittarius (Sgr) stre am assuredly prohibit str ongly flattened dark 
haloes with q < 0.6 dEvans & Bowdenll2014l) . Whilst a definitive 
answer from Sgr stream kinematics must await a thorough explo¬ 
ration of stream generation in flattened haloes using modem algo¬ 
rithms dGibbons et alll201^ . the debris of the Sgr stretches a full 
360° over the sky and is almost confined to a plane. 

The self-consistent distribution functions (DFs) in terms of 
the energy E and angular momentum L have been given for the 
models with constant anisotropy. There are however strong reasons 
for usin g actions, inst ead of integrals of t he motion like energy , 
in DFs dBinnevll2013h . Here, we note that IWilliams et all d201'^ 
have shown that the Hamiltonian as a function of actions in scale- 
free power-laws is an almost linear function of the actions, en¬ 
abling schemes to be easily devised to c onver t the DFs into ac¬ 
tion space if desired. lEvans & William^ d2014b provide a practi¬ 
cal exa mple of such an a lgori thm for one member of th e doubloon 
family. fPosti et al.ld2015h and I Williams & Evan3 d2015h give algo¬ 
rithms for action-based DFs for models with density falling like 
power-laws in certain regimes. 


The self-consistent DFs describe the velocity distributions of 
dark matter needed to sustain the doubloon models. For applica¬ 
tions to stars in stellar halos or in elliptical galaxies, we are in¬ 
terested in the properties of luminous tracers within the doubloon 
models. We have provided a number of results for both power-law 
and Einasto tracer populations. There is good evidence that the ve¬ 
locity dispersion tensor of th e Milky Way’s stellar halo is aligned i n 
spherical polar coordinates dSmith et al]|20093 : [Bond et al.l2010l) . 
Additionally, the velocity ellipsoid of stellar populations in ellipti¬ 
cal galaxies is probably aligned in spheroidal coordinates, which 
asymptotically become spherical. Therefore, spherically aligned 
Jeans solutions are of considerable astrophysical interest. We iden¬ 
tify a simple algorithm for solving the spherically aligned Jeans 


equations for flattened tracers in spherical potentials, and pro¬ 
vide solutions for power-law and hyper-Einasto tracers (defined in 
eq.l^. 

The distribution of high velocity stars of the stellar halo of the 
Milky Way are becoming available dSmith et al.ll2007l : IPifS et ^ 
|20l4lHawkins et al.ll2015h . So, we have provided the asymptotic 
forms of the DFs for tracer populations with power-law and Einasto 
density distributions. The form of the high velocity tail is particu¬ 
larly interesting as this may betray properties of the dark matter 
halo. Power-law tracers have velocity distributions with power-law 
or Maxwellian tails. Einasto tracers have super-exponential cut-offs 
in the velocity distributions. Although the observational data are of¬ 
ten fitted to models like f(v) oc - vY this is only strictly correct 
for power-law tracers. If the stellar halo is described by an Einasto 
profile, then 


(fiv) oc (jJesc - vf exp[- ^^ , ^ 

[/(u) = 0 V> Uesc, 


(49) 


with A, B and C constants, is a better description of the high veloc¬ 
ity tail near the escape speed. In principle, the different forms of the 
velocity distributions of tracers can provide us with evidence on the 
extent of the dark matter potential. This is a subject on which there 
is not merely little hard evidence, but very few avenues in which to 
gain evidence. 
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APPENDIX A: CONSTANT ANISOTROPY MODELS 
AI F ractional calculus 


Many of the formulae concerning the constant anisotropy DR can 
be swiftly derived using the notion of fraction al differentiation. Le t 
us consider the Riemann-Lioville integral fsee lErdelvi et al.ll954h : 

:TJ(x) (M > 0). (Al) 

This generalizes Cauchy’s formula for repeated integration for an 
arbitrary positive order. This is also recognized as (generalized) 
Abel transform for 0 < yu < 1 with the classical case resulting 
from /J = 1/2. These obey the composition rule 

= (A2) 


and the differentiation follows 



:irm 

fix) 


in < ^) 

in =u) 


(A3) 


Thus the Riemann-Lioville integral may be inverted through 
gix) = :fj{x) => fix) = (-) llXgix) (A4) 


where [juJ and e = yU - [juJ are the integer floor and the fraction part 
of /i, respectively. Let us also define the fractional derivative: 


a^Jfix) = 


r(^+i) 

27ti 


(-V+) 


I 


fiz) dz 
(z - x)'+f 


a 


(A5) 


by means of the complex contour integral. Here, the contour starts 
and ends at the base point a, and encircles x in the counter¬ 
clockwise direction. Thanks to Cauchy’s integral (and differenti¬ 
ation) formula, this coincides with the customary notion of dif¬ 
ferentiation for integer values of that is, fix) = fix) and 
a^x'fix) = T"\x) for n e {1,2,...). Explicit calculations can 
demonstrate that the fractional derivative is the inverse operator of 
the Riemann-Lioville integral, whereby 

.dXjjix) = fix), (A6) 

and so follows that a^x^fix) = iA!Axf^XP ^fi^) where \f~\ is the 
integer ceiling of f and also assuming ^f^fix) - fix). 


A2 Constant anisotropy DF 

Suppose that the DF is given by the ansatz of equation (dill. The 
density profile results from the integral over the velocity, 

FiE, D) dEdD 


P = 




E>a,O>0,v)>0 XfiEjJxy) 


(A7) 


where and a = firf) if ro is the finite boundary radius or a = 
limr_>oo fir). Hence if we define the augmented density 


gif) ^ ildfp = llj-^Mf) 


(A8) 
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the energy part of the DF is inverted as 

3 / d 

/£(£) = = W"ci{E), (A9) 

where n = [(3/2-//) | and e = 'ill —ji — n , which may be compared 
with equation (3) of lEvans & AnI bOOfili (note the scale constants 
are chosen differently here). For the isotropic case (yS = 0), we have 

which reproduces Eddington’s (1916) formula, which can be 
thought of as the fractional derivative of order 3/2. 

Eor the DE of the form of equation d, the energy part /fifE) 
is directly related to the local energy distribution 

dp fEiE) 

dE F(3/2-;S) (2r^f ^ ’ 

as well as the speed distribution 
dp v^-VfE[,lj{r)-v^l2} 

dt)2 23/2F(3/2-/?)r2f* ■ ^ ’ 

Flere also note that the local escape speed is given by - yj2ij/{r) 
provided that limr^^oo i/'(r) = 0. 


and so equation l lA9b results in the DE in the form of equation ( 112b 
given by equation G3}. Eor models on the inner branch, we have 

C(l-ip|e)(l-e)f+‘ 


= g(s ), g(s) = 


(A16) 


where ^ = 2\p\ '(1 - /I) > 0. If// = | - n, equation l lA9b becomes 


fE{E) = 


d"g(E) 

dE" 




(A17) 


e=exp(-|;)|£/V2) 

which can be computed analytically. For others, we can still expand 
g{s) in a power series in e (note 0 < s < 1 since p < 0 but ij/ > Qi)\ 


gW = 


C(1 - \p\s) 


z 

1=0 


'^+1 

j 


{-sy. 


(A18) 


Note, if ^ is a non-negative integer, the sum terminates after the 
7 = ^ -I- 1 term and so reduces to a polynomial in e. Given equation 
( IA14b . applying equation dA9b then results in /^(E) given as a sum 
over incomplete gamma functions, P\f} - |,-(y - ^)|p|£] where 
£ = E/V^(> 0). If p is an integer, these are reducible to the error 
function or Dawson’s integral, particular examples of which are 
provided in the main body. 


A3 Auxiliary results 

For models on the outer branch, we make us e of Hankel’s Loop 
integral f or the reciprocal g amma function (see lErdelvi et al.ll953L 
ea.l.6(2): l01ver et al ]|20ia § 5.9(1)) indicating 

(E+) 

For models on the inner branch however, we need to intro¬ 
duce the incomplete gamma function. The contour i ntegral repre¬ 
senta tion of the lower incomplete g amma function (cf. lErdelvi et al.l 
Il95i eQ.9.3(l): loiver et al.ll201(l § 8.6(2)) implies 

(£+) 

g ^ nd + 1) r exp(.?i/r)di/' 

“ ® “ 27ti J (iA-E)‘+'i (A14) 

0 

= i''e'-^E(-T, sE) = E-''e*^y*(-T, sE) 

Here, P(a,x) = y(a, x)IT(a) = x‘‘y'‘(a,x) is the regularized 
lower incomplete gamma function. If T is a positive integer, then 
P(-A, x) = I and so this is same as the ordinary n-th derivative. 
For a half-integer A, this reduces to the error function (or Dawson’s 
integral, which is equivalent to the imaginary error integral) - note 
erf(x) = (1/ V^)7(1/2 ,x^) = E(l/2,x^) anderfi(x) = -ierf(ix). 


A4 Constant anisotropy DFs of self-consistent Doubloons 


For models on the outer branch, p > 0 and (ft < 0 and so 0 < 
e < 1 where e = exp(pif//V^). Then of equation dl lb may be 
expanded to the power-series of e using binomial series; 


gW = 


C(e + pe^) ^ (Ay 

g-r Zj ;i ^ 

j=o J- 


= C 




r- 


(y-1)!/ 


(A15) 


where A = 2p *(1 -/!) - 1 and C = 2^ KnG). Equation 

dA13b then indicates = V^^^[p(A+ 1 -l-y)] 
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